##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#larv.ranef <- read.csv("larv.mod.ranef.csv")
larv.ranef <- read.csv("larv.mod.ranef.all.csv")
##ps object & taxa table
#ps.lar.trim <- readRDS("ps.exp.lar.trim.rds")
ps.lar.trim <- readRDS("ps.exp.all.lar.trim.rds")
tax.lar.trim <- data.frame(ps.lar.trim@tax_table)## Loading required package: phyloseq
##individual ASVs in all treats
##raw counts
ps.lar.trim.asv <- subset_taxa(ps.lar.trim,id=="ASV0001")
plot_bar(ps.lar.trim.asv)+
facet_wrap(Food_type~Microbe_treatment,scales="free",ncol=5)ps.lar.trim.rel = transform_sample_counts(ps.lar.trim, function(x) x / sum(x))
#ps.lar.trim.rel.trt <- subset_samples(ps.lar.trim.rel,Microbe_treatment=="MEM")
ps.lar.trim.rel.asv <- subset_taxa(ps.lar.trim.rel,id=="ASV0001")
plot_bar(ps.lar.trim.rel.asv)+
facet_wrap(Food_type~Microbe_treatment,ncol=5,scales="free")##
## variable variable:Food_type
## 54 108
## variable:Microbe_treatment variable:Microbe_treatment:Food_type
## 270 540
#table(larv.ranef$term)
larv.ranef.treat <- subset(larv.ranef,grpvar=="variable:Microbe_treatment:Food_type")
larv.ranef.treat$id <- substr(larv.ranef.treat$grp,1,7)
larv.ranef.treat$food <- str_sub(larv.ranef.treat$grp,13,18)
larv.ranef.treat$microbes <- str_sub(larv.ranef.treat$grp,9,11)
larv.ranef.treat$food_micr <- paste0(larv.ranef.treat$food,"_",larv.ranef.treat$microbes)
#larv.ranef.treat$food_micr <- substr(larv.ranef.treat$grp,9,18)
##just the columns for stuff below:
larv.ranef.less <- larv.ranef.treat %>%
select(-X,-term,-grp,-condsd,-food,-microbes)
# Reshape the data frame to have the food types as two columns
larv.ranef.wide <- larv.ranef.less %>%
pivot_wider(
names_from = food_micr, # Column to use for new column names
values_from = condval # Column to fill the new columns
)
##add taxa info
larv.ranef.wide.tax <- merge(larv.ranef.wide,tax.lar.trim,by="id")
##more informative ASV ids
larv.ranef.wide.tax$asv.gen <- paste0(larv.ranef.wide.tax$id,"_",larv.ranef.wide.tax$Genus)library("pheatmap")
library("RColorBrewer")
cond.df <- larv.ranef.wide.tax[,c("asv.gen","Bamboo_EMO","Bamboo_MMO","Bamboo_MEM","Larval_EMO","Larval_MMO","Larval_MEM")]
#subset(cond.df)
row.names(cond.df) <- cond.df$asv.gen
cond.df2 <- cond.df[,2:7]
cond.mat <- as.matrix(cond.df2)
#row.names(otu.rel) == row.names(sam.rel)
pheatmap(cond.mat,
#cluster_cols=F,
#color = colorRampPalette(c("blue", "white", "red"))(100),
color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)pheatmap(cond.mat,
scale="row",
#cluster_cols=F,
#color = colorRampPalette(c("blue", "white", "red"))(100),
color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)##MMO - MMO BETWEEN FOODS
larv.ranef.wide.tax$diff.bmmo.lmmo <- larv.ranef.wide.tax$Bamboo_MMO - larv.ranef.wide.tax$Larval_MMO
##EMO - EMO BETWEEN FOODS
larv.ranef.wide.tax$diff.bemo.lemo <- larv.ranef.wide.tax$Bamboo_EMO - larv.ranef.wide.tax$Larval_EMO
##MEM - MEM BETWEEN FOODS
larv.ranef.wide.tax$diff.bmem.lmem <- larv.ranef.wide.tax$Bamboo_MEM - larv.ranef.wide.tax$Larval_MEM
##All & eco
##EMO - EMO BETWEEN FOODS
larv.ranef.wide.tax$diff.ball.lall <- larv.ranef.wide.tax$Bamboo_ALL - larv.ranef.wide.tax$Larval_ALL
##MEM - MEM BETWEEN FOODS
larv.ranef.wide.tax$diff.beco.leco <- larv.ranef.wide.tax$Bamboo_ECO - larv.ranef.wide.tax$Larval_ECO
##MMO - EMO BETWEEN FOODS
# larv.ranef.wide.tax$diff.bmmo.lemo <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Bamboo_EMO
# ##MMO - MEM BETWEEN FOODS
# larv.ranef.wide.tax$diff.bmmo.lmem <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Bamboo_MEM
##EMO - MMO BETWEEN FOODS
#larv.ranef.wide.tax$diff.lemo.bmmo <- larv.ranef.wide.tax$Larval_EMO - larv.ranef.wide.tax$Bamboo_MMO
##EMO - MEM BETWEEN FOODS
#larv.ranef.wide.tax$diff.lemo.bmem <- larv.ranef.wide.tax$Larval_EMO - larv.ranef.wide.tax$Bamboo_MEM
##MEM - MMO BETWEEN FOODS
#larv.ranef.wide.tax$diff.lmem.bmmo <- larv.ranef.wide.tax$Larval_MEM - larv.ranef.wide.tax$Bamboo_MMO
##MEM - EMO BETWEEN FOODS
#larv.ranef.wide.tax$diff.lmem.bemo <- larv.ranef.wide.tax$Larval_MEM - larv.ranef.wide.tax$Bamboo_EMO##MMO - EMO WITHIN BAMBOO
#larv.ranef.wide.tax$diff.bmmo.bemo <- larv.ranef.wide.tax$Bamboo_MMO - larv.ranef.wide.tax$Bamboo_EMO
##MMO - MEM WITHIN BAMBOO
larv.ranef.wide.tax$diff.bmmo.bmem <- larv.ranef.wide.tax$Bamboo_MMO - larv.ranef.wide.tax$Bamboo_MEM
##EMO - MEM WITHIN BAMBOO
larv.ranef.wide.tax$diff.bemo.bmem <- larv.ranef.wide.tax$Bamboo_EMO - larv.ranef.wide.tax$Bamboo_MEM
##MMO - EMO WITHIN LARVAL FOOD
#larv.ranef.wide.tax$diff.lmmo.lemo <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Larval_EMO
##MMO - MEM WITHIN LARVAL FOOD
larv.ranef.wide.tax$diff.lmmo.lmem <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Larval_MEM
##EMO - MEM WITHIN LARVAL FOOD
larv.ranef.wide.tax$diff.lemo.lmem <- larv.ranef.wide.tax$Larval_EMO - larv.ranef.wide.tax$Larval_MEMlarv.ranef.diffs <- larv.ranef.wide.tax %>%
select(id, starts_with("diff.")) %>%
pivot_longer(
cols = starts_with("diff."),
names_to = "comparison",
values_to = "difference"
)
#larv.ranef.diffs.asv1 <- subset(larv.ranef.diffs,id=="ASV0001")
#ggplot(larv.ranef.diffs.asv1,aes(x=comparison,y=difference))+
# geom_bar(stat="identity")+
# theme(axis.text.x=element_text(angle=45,hjust=1))
##alternatively
ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bemo.lemo),y=diff.bemo.lemo))+
geom_bar(stat="identity")+
theme(axis.text.x=element_text(angle=45,hjust=1))+
ggtitle("EMO - bamboo minus larval food")##individual ASVs within treatments
ps.lar.trim.rel = transform_sample_counts(ps.lar.trim, function(x) x / sum(x))
ps.lar.trim.rel.trt <- subset_samples(ps.lar.trim.rel,Microbe_treatment=="MEM")
ps.lar.trim.rel.trt.asv <- subset_taxa(ps.lar.trim.rel.trt,id=="ASV0027")
plot_bar(ps.lar.trim.rel.trt.asv)+
facet_wrap(Food_type~Microbe_treatment)ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmmo.lmmo),y=diff.bmmo.lmmo))+
#geom_bar(stat="identity",position=position_dodge(1))+
geom_bar(stat="identity")+
#scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
ggtitle("MMO - bamboo minus larval food")ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bemo.lemo),y=diff.bemo.lemo))+
#geom_bar(stat="identity",position=position_dodge(1))+
geom_bar(stat="identity")+
#scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
ggtitle("EMO - bamboo minus larval food")ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmem.lmem),y=diff.bmem.lmem))+
#geom_bar(stat="identity",position=position_dodge(1))+
geom_bar(stat="identity")+
#scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
ggtitle("MEM - bamboo minus larval food")library("pheatmap")
library("RColorBrewer")
#diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.beco.leco","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem","diff.ball.lall")]
diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem")]
row.names(diff.df) <- diff.df$asv.gen
diff.df2 <- diff.df[,2:4]
#diff.df2 <- diff.df[,2:6]
diff.mat <- as.matrix(diff.df2)
#row.names(otu.rel) == row.names(sam.rel)
pheatmap(diff.mat,
color = colorRampPalette(c("#8C510A", "white", "darkolivegreen4"))(100),
scale="row",
cluster_rows=T,
cluster_cols=T,
#color = colorRampPalette(brewer.pal(n=7,name="BrBG"))(100)
#color = colorRampPalette(brewer.pal(n=7,name="PiYG"))(100)
#color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)## [1] "#8C510A"
ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmmo.lmmo),y=diff.bmmo.lmmo))+
#geom_bar(stat="identity",position=position_dodge(1))+
geom_bar(stat="identity")+
#scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
ggtitle("MMO - bamboo minus larval food")ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bemo.lemo),y=diff.bemo.lemo))+
#geom_bar(stat="identity",position=position_dodge(1))+
geom_bar(stat="identity")+
#scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
ggtitle("EMO - bamboo minus larval food")ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmem.lmem),y=diff.bmem.lmem))+
#geom_bar(stat="identity",position=position_dodge(1))+
geom_bar(stat="identity")+
#scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
theme_cowplot()+
theme(axis.text.x=element_text(angle=90,hjust=1))+
ggtitle("MEM - bamboo minus larval food")library("pheatmap")
library("RColorBrewer")
diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.beco.leco","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem","diff.ball.lall")]
#diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem")]
row.names(diff.df) <- diff.df$asv.gen
#diff.df2 <- diff.df[,2:4]
diff.df2 <- diff.df[,2:6]
diff.mat <- as.matrix(diff.df2)
#row.names(otu.rel) == row.names(sam.rel)
pheatmap(diff.mat,
color = colorRampPalette(c("#8C510A", "white", "darkolivegreen4"))(100),
scale="row",
cluster_rows=T,
cluster_cols=T,
#color = colorRampPalette(brewer.pal(n=7,name="BrBG"))(100)
#color = colorRampPalette(brewer.pal(n=7,name="PiYG"))(100)
#color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)## [1] "#8C510A"
#install.packages("ggeffects")
library("ggeffects")
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models")
larv.mod.nb1 <- readRDS("larv.mod.nb1.rds")
#ps.trim.wtr <- readRDS("ps.exp.wtr.trim.rds")
#ps.trim.wtr
ps.trim.lar <- readRDS("ps.exp.lar.trim.rds")
ps.trim.larpred.lar.nb1 <- ggpredict(larv.mod.nb1, terms=c("variable", "Food_type","Microbe_treatment"), type="re")
plot(pred.lar.nb1)+
theme(axis.text.x=element_text(angle=90))+
theme_classic()
# pred.lar.nb2 <- ggpredict(larv.mod.nb2, terms=c("variable", "Food_type","Microbe_treatment"), type="re")
# plot(pred.lar.nb2)
tax.trim.lar <- data.frame(ps.trim.lar@tax_table)
##make same column names
pred.lar.nb1$id <- pred.lar.nb1$x
pred.lar.tax <- merge(tax.trim.lar,pred.lar.nb1,by="id")
ggplot(pred.lar.tax,aes(x=group,y=predicted,color=facet))+
geom_point()+
facet_wrap(~facet)
pred.lar.tax.mmo <- subset(pred.lar.tax,facet=="MMO")
ggplot(pred.lar.tax.mmo,aes(x=group,y=predicted,fill=id))+
#geom_bar(stat="identity",position=position_dodge(0.5))+
geom_bar(stat="identity")+
facet_wrap(~Genus,scales="free")+
theme(legend.position="none")
pred.lar.tax.bamb <- subset(pred.lar.tax,group=="Bamboo")
ggplot(pred.lar.tax.bamb,aes(x=facet,y=predicted,fill=id))+
geom_bar(stat="identity")+
facet_wrap(~Family,scales="free")+
theme(legend.position="none")Trying some different plots
library("cowplot")
pred.lar.tax.comp <- subset(pred.lar.tax,facet!="ECO"&facet!="ALL")
pred.lar.tax.comp.bam <- subset(pred.lar.tax.comp,group=="Bamboo")
ggplot(pred.lar.tax.comp.bam,aes(x=facet,y=predicted,fill=facet))+
facet_wrap(~id,scales="free")+
geom_bar(stat="identity")+
theme_cowplot()
##heatmap
library("pheatmap")
library("RColorBrewer")
library("tidyr")
#pred.lar.tax.comp$asv_food <- paste0(pred.lar.tax.comp$id,"_",pred.lar.tax.comp$group)
colnames(pred.lar.tax.comp)
##remove not needed columns
pred.lar.min <- pred.lar.tax.comp %>%
select(predicted,facet,group,id)
##bamboo
pred.lar.min.bam <- subset(pred.lar.min,group=="Bamboo")
pred.lar.min.bamw <- pred.lar.min.bam %>%
pivot_wider(names_from=id,values_from=predicted) %>%
as.data.frame()
row.names(pred.lar.min.bamw) <- pred.lar.min.bamw$facet
pred.lar.min.bamw2 <- pred.lar.min.bamw[,3:ncol(pred.lar.min.bamw)]
pheatmap(as.matrix(pred.lar.min.bamw2),
color = colorRampPalette(c("white", "orange", "red"))(100),
#color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100),
#scale="column",
cluster_cols = F,
cluster_rows=F
)
##larval food
pred.lar.min.lar <- subset(pred.lar.min,group=="Larval food")
pred.lar.min.larw <- pred.lar.min.lar %>%
pivot_wider(names_from=id,values_from=predicted) %>%
as.data.frame()
row.names(pred.lar.min.larw) <- pred.lar.min.larw$facet
pred.lar.min.larw2 <- pred.lar.min.larw[,3:ncol(pred.lar.min.larw)]
pheatmap(as.matrix(pred.lar.min.larw2),
color = colorRampPalette(c("white", "orange", "red"))(100),
#color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100),
#scale="column",
cluster_cols = F,
cluster_rows=F
)
library("ggrepel")
##correlating what happened in bamboo & larval food
pred.lar.min$group <- gsub("Larval food","Larval",pred.lar.min$group)
pred.lar.min.cor <- pred.lar.min %>%
pivot_wider(names_from=group,values_from=predicted)
pred.lar.min.cor$facet <- factor(pred.lar.min.cor$facet,levels=c("MMO","EMO","MEM"))
##full plot
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
geom_point(size=2)+
geom_abline(intercept=0,slope=1,linetype="dotted")+
geom_text_repel(color="black")+
theme_cowplot()+
ylab("Standard diet")+
xlab("Bamboo leaves")+
scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))
##zooming into the low section
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
geom_point()+
geom_abline(intercept=0,slope=1,linetype="dotted")+
geom_abline(intercept=0,slope=0.5,linetype="dotted")+
geom_abline(intercept=0,slope=2,linetype="dotted")+
xlim(0,100)+
ylim(0,100)+
geom_text_repel(color="black")+
theme_cowplot()+
ylab("Standard diet")+
xlab("Bamboo leaves")+
scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))
mean(sample_sums(ps.lar.trim))
#47762.4
1000/47762.4
##max 2k - larval food
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
geom_point()+
geom_abline(intercept=0,slope=1,linetype="dotted")+
geom_abline(intercept=0,slope=0.5,linetype="dotted")+
geom_abline(intercept=0,slope=2,linetype="dotted")+
#xlim(0,4000)+
ylim(1000,30000)+
geom_text_repel(color="black")+
theme_cowplot()+
ylab("Standard diet")+
xlab("Bamboo leaves")+
scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))
##max 2k - bamboo
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
geom_point()+
geom_abline(intercept=0,slope=1,linetype="dotted")+
geom_abline(intercept=0,slope=0.5,linetype="dotted")+
geom_abline(intercept=0,slope=2,linetype="dotted")+
xlim(1000,30000)+
#ylim(1000,30000)+
geom_text_repel(color="black")+
theme_cowplot()+
ylab("Standard diet")+
xlab("Bamboo leaves")+
scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))Checking what the actual data look like
ps.trim.lar.comp <- subset_samples(ps.lar.trim,Microbe_treatment!="ALL"&Microbe_treatment!="ECO")
ps.trim.lar.comp
ps.lar.trim.rel = transform_sample_counts(ps.trim.lar.comp, function(x) x / sum(x))
ps.lar.trim.rel
ps.asv1 <- subset_taxa(ps.trim.lar.comp,id=="ASV0001")
ps.asv1@otu_table
plot_bar(ps.asv1,x="Sample")+
facet_wrap(Microbe_treatment~Food_type)
plot_bar(ps.asv1)+
facet_wrap(Microbe_treatment~Food_type,scales="free")
ps.asv3 <- subset_taxa(ps.lar.trim.rel,id=="ASV0003")
ps.asv3
plot_bar(ps.asv3,x="Microbe_treatment")+
facet_wrap(~Food_type)
plot_bar(ps.asv3)+
facet_wrap(Microbe_treatment~Food_type,scales="free")
ps.asv4 <- subset_taxa(ps.lar.trim.rel,id=="ASV0004")
ps.asv4
plot_bar(ps.asv4,x="Microbe_treatment")+
facet_wrap(~Food_type)
ps.asv6 <- subset_taxa(ps.lar.trim.rel,id=="ASV0006")
ps.asv6
plot_bar(ps.asv6,x="Microbe_treatment")+
facet_wrap(~Food_type)
ps.asv7 <- subset_taxa(ps.lar.trim.rel,id=="ASV0007")
ps.asv7
plot_bar(ps.asv7,x="Microbe_treatment")+
facet_wrap(~Food_type)
ps.asv19 <- subset_taxa(ps.trim.lar.comp,id=="ASV0019")
ps.asv19
plot_bar(ps.asv19,x="Microbe_treatment")+
facet_wrap(~Food_type)
plot_bar(ps.asv19)+
facet_wrap(Microbe_treatment~Food_type,scales="free")
ps.asv33 <- subset_taxa(ps.trim.lar.comp,id=="ASV0033")
ps.asv33
plot_bar(ps.asv33,x="Microbe_treatment")+
facet_wrap(~Food_type)
plot_bar(ps.asv33)+
facet_wrap(Microbe_treatment~Food_type,scales="free")
ps.asv54 <- subset_taxa(ps.lar.trim.rel,id=="ASV0054")
ps.asv54
plot_bar(ps.asv54,x="Microbe_treatment")+
facet_wrap(~Food_type)
plot_bar(ps.asv54)+
facet_wrap(Microbe_treatment~Food_type)
larv.ranef.wide.tax[6,]
# )+
# facet_wrap(Food_type~Microbe_treatment)